Gauge conditions for long-term numerical black hole evolutions without excision 
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Numerical relativity has faced the problem that standard 3+1 simulations of black hole spacetimes 
without singularity excision and with singularity avoiding lapse and vanishing shift fail after an 
evolution time of around 30 — 40M due to the so-called slice stretching. We discuss lapse and shift 
conditions for the non-excision case that effectively cure slice stretching and allow run times of 
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I. INTRODUCTION 

A crucial role in numerical relativity simulations of 
two black holes (BH) is played by the choice of coor- 
dinates. This gauge choice involves both the choice of 
a lapse function and of a shift vector, which typically 
have to be determined dynamically during a numerical 
evolution. The first results for colliding BH's were ob- 
tained for head-on collisions using the ADM decomposi- 
tion of the Einstein equations with the lapse determined 
by the maximal slicing condition and the shift vector set 
to zero • Maximal slices are known to be singular- 
ity avoiding, that is, starting from BH initial data where 
the physical singularity is to the future of the initial hy- 
persurface, the lapse approaches the Minkowski value of 
unity in the asymptotically flat regions, but approaches 
zero near the physical singularity. In this way one can 
in principle foliate a BH spacetime without singularities, 
but since time marches on in the far regions while being 
frozen in the interior, the slices become more and more 
distorted. Historically, this phenomenon has been called 
'grid stretching' by the numerical relativity community, 
though we will refer to it as 'slice stretching' since it is 
a property of the slices themselves, quite independent of 
the existence of a numerical grid. Slice stretching intro- 
duces a difficult problem for numerical simulations since 
the metric develops large gradients that keep on growing 
until the numerical code can no longer handle them and 
fails. Advanced numerical methods can help in spherical 
symmetry, see, e.g., Q, but to date have not proved very 
successful in three-dimensional (3D) evolutions ||. 

Nonetheless, such singularity avoiding slicings with 
vanishing shift have been quite successful, since they do 
allow a finite evolution time of roughly t = 30 — 40M 
(with M the total ADM mass of the system) for fully 
3D evolutions of BH's, as first demonstrated in 1995 for 
the case of a sing le Schwarzschild BH §]. In @ the first 
fully 3D simulation of the grazing collision of two nearby 
BH's was performed with singularity avoiding slicing and 



vanishing shift, lasting for about 7M. With improved 
techniques the grazing collision has recently been pushed 
to about 35M, which for the first time allowed the ex- 
traction of gravitational waveforms from a 3D numerical 
merger ||. And even though singularity avoiding slic- 
ings with vanishing shift have so far been limited to a 
finite time interval of only 30 — 40M, this interval can 
be moved into the truly non-linear regime of a plunge 
starting from an approximate innermost stable circular 
orbit (ISCO) of two BH's, since the remainder of the 
merger and ring-down can be computed using the close 
limit approximation || . Following such an approach, the 
first waveforms for the plunge from an approximate ISCO 
have been obtained ]Tc|-[T2^| . 

So far the most important strategy to avoid slice 
stretching has been black hole excision |l^JT^]. The idea 
is to use horizon penetrating coordinates (notice that 
maximal slicing is horizon penetrating unless one imposes 
the extra boundary condition of having a vanishing lapse 
at the horizon) and to excise the interior of the BH's from 
the numerical grid. A non-vanishing shift is essential to 
keep grid points from falling into the BH. This approach 
has seen many successful implementations for single black 
holes. First demonstrated in 3D in ||, with further de- 
velopment it has, in particular, allowed to move a black 
hole across the numerical grid . If a stable numerical 
implementation can be found, this approach should make 
it possible to simulate many orbits of two well separated 
BH's. The key difference between BH excision and the 
use of singularity avoiding slicings with vanishing shift is 
that with excision single static BH's can be stably evolved 
for thousands of M; see for the case of evolutions us- 
ing null coordinates (that do not directly generalize to 
binary BH systems), and recently 100, 000A7 have been 
reached for a single BH with a 3+1 Cauchy code in oc- 
tant symmetry [[16]]. Black hole excision holds a lot of 
promise, even if currently evolutions of only 9-15M have 
been achieved for binary BH's [ p7| . 

In this paper we demonstrate that the new lapse and 



1 



shift conditions introduced in for the case of a single 
distorted BH with excision (using the excision techniques 
of ]l6| ) can work well even without excision. This allows 
us to break through the barrier of about 35Af for singu- 
larity avoiding slicings in 3+1 numerical relativity. Our 
gauge choice maintains singularity avoidance but cures 
the main problems associated with slice stretching, al- 
lowing us to reach 500M and more for the evolution of 
single or even distorted BH's. For BH's colliding head- 
on that merge early on during the evolution, i.e. which 
start out sufficiently close to each other, the final BH can 
again be evolved for hundreds or even thousands of M. 

Moreover, these gauge conditions have two impor- 
tant effects: (a) they drive the system towards a static 
state, virtually, if not completely, eliminating the chronic 
growth in metric functions typical of slice stretching. 
Hence, in principle they should allow for indefinitely long 
evolutions (if no other instabilities develop; see below), 
(b) since unbounded growth in metric functions is halted, 
they allow much more accurate results to be obtained for 
extremely long times, and at lower resolution than be- 
fore. Below we will show results obtained for colliding 
black holes that show only 10% error in the horizon mass 
after more than 5000M of evolution. 

The evolutions in this paper are carried out using the 
puncture method for evolutions [0J^| . The starting point 
is initial data of the Brill-Lindquist topology in isotropic 
coordinates JlSfl - This 'puncture' data is defined on R 3 
minus a point for each of the internal asymptotically flat 
ends of the BH's. If one treats the coordinate singular- 
ity at the punctures appropriately, the punctures do not 
evolve as long as the shift vanishes there. That is, the 
metric and the extrinsic curvature do not evolve at the 
punctures. It can also be checked that the maximal slic- 
ing equation produces a smooth numerical solution for 
the lapse at the punctures. 

One basic observation for our choice of shift vector is 
that the "Gamma freezing" shift introduced in |l^] for 
our project in simple BH excision has the following prop- 
erty when the BH's are not excised but are represented 
as punctures: Initially the shift is zero, but as the slice 
stretching develops, the shift reacts by pulling out points 
from the inner asymptotically flat region near the punc- 
tures. The lapse and shift conditions taken together are 
then able to virtually stop the evolution of one or even 
two black holes, essentially mimicking the behavior of the 
lapse and shift known from stable evolutions of a BH in 
Kerr-Schild coordinates. This is a key result that will be 
detailed below. 

The paper is organized as follows. First we introduce 
the evol utio n equations and the constraints in Section ||. 
In Sec. |n| we discuss the gauge conditions. The punc- 
ture initial data and puncture evolutions are discussed in 
Sec. |y| and Sec. [v|. In Sec. VI miscellaneous aspect s of 
our numerical implementation are discussed. In Sec. VII 
we present numerical results for one and two BH's, and 
we conclude in Sec. VIII. 



II. FORMULATION 

The standard variables in the 3+1 formulation of ADM 
(Arnowitt-Deser-Misner, see [Q) are the 3- metric 7y 
and the extrinsic curvature K^ . The gauge is determined 
by the lapse function a and the shift vector (3 l . We will 
only consider the vacuum case. The evolution equations 
are 



(d t - Cp) K i:j 



D.Dja 



+a(Rij 
and the constraints are 



KKi. 



2K ik K k j), 



H = R 

V 1 = 



K - K tj K l] = 0, 



7 ij K) = 0. 



(1) 
(2) 



(3) 
(4) 



Here Cp is the Lie derivative with respect to the shift vec- 
tor (3 1 , Di is the covariant derivative associated with the 
3-metric 7^ , Rij is the three-dimensional Ricci tensor, R 
the Ricci scalar, and K is the trace of Kij. 

We will use the BSSN form of these equations (Baum- 
garte, Shapiro p]]] , and Shibata, Nakamura Q). One 
introduces new variables based on a trace decomposition 
of the extrinsic curvature and a conformal rescaling of 
both the metric and the extrinsic curvature. The trace- 
free part of the extrinsic curvature is defined by 



Aij = - 



(5) 



Assuming that the metric 7^ is obtained from a confor- 
mal metric jij by a conformal transformation, 



lij = "0 4 7y> 



(0) 



we can choose a conformal factor ip such that the deter- 
minant of "fij is 1: 



V-7 1/12 , 

lij V 7 lij 7 ^ lij 1 

7=1, 



(7) 

(8) 
(9) 



where 7 is the determinant of 7^ and 7 is the determinant 



of 



lij- 



Instead of 7^ and Ku we can therefore use the 



variables 



= = ^2 ln 7, 



K 

% 
Ai 3 



lij 



K lJ , 



e^Hj, 



-40 A 



13 1 



(10) 

(11) 
(12) 

(13) 



where 7^ has determinant 1 and A^ has vanishing 
trace. Furthermore, we introduce the conformal connec- 
tion functions 
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r = ^ k r jk = -d^, (w) 

where T*j k is the Christoffel symbol of the conformal 
metric and where the second equality holds only if the 
determinant of the conformal 3-metric 7 is unity (which 
is true analytically but may not hold numerically). We 
call </>, K, 7y, Aij, and T l the BSSN variables. 

In terms of the BSSN variables the evolution equation 
(|l|) becomes 

{dt-Cp)% =-2aAij, (15) 
(8t-C fi )4, = ~aK, (16) 

while equation (||) leads to 

{d t - Cp) Aij = e-^l-DiDja + aRij] TF 

+ a(KA l3 - 2A ik A k j), (17) 
[dt -Cp)K= -D l D t a 

+ a(A ii & + \k 2 ), (18) 

where TF denotes the trace-free part of the expression 
in brackets with respect to 7^. Note that the right- 
hand side of the evolution equation ( |17|) for the trace- 
free variable Aij is trace-free except for the term Ai k A k j. 
This is no contradiction since the condition that A^ 
remains trace-free is (dt — Cp) (^ Aij) = and not 
fi{d t -Lp)Aij = Q. 

On the right-hand side of (|lj) we have used the Hamil- 
tonian constraint (||) to eliminate the Ricci scalar, 

R = A,,/\ ■' - K 2 = A. l:j A zj - ^K 2 . (19) 

The momentum constraint ^ becomes 

djA ij = -f l 3k A 3k - <oA ij dj<t> + ~j ij djK. (20) 

An evolution equation for T l can be obtained from 
and @, 

d t f l = -2 (adjA ij + A»dja\ - dj (C f j ) , (21) 

where we will use the momentum constraint above to 
substitute for the divergence of A l K One subtlety in 
obtaining numerically stable evolutions with the BSSN 
variables is precisely the question of how the constraints 
are used in the evolution equations. Several choices are 
possible and have been studied, see e.g. ]23| ]. 

Note that in the preceding equations we are computing 
Lie derivatives of tensor densities. If the weight of a 
tensor density T is w, i.e. if T is a tensor times r ) w ^ 2 , 
then 

C T=[C T\T°+wT8 k k , (22) 



where the first term denotes the tensor formula for Lie 
derivatives with the derivative operator d and the second 
is the additional contribution due to the density factor. 
The density weight of if; — is i , so the weight of 7^ 

and A^ is — I and the weight of 7^ is |. To be explicit, 

C <j> = k 8 k <f>+±8 k k , (23) 
Cp% = k 8 k % + % k 8j0 k + j jk di0 k - \lijd k k , (24) 
C/3f j = k 8 k f j - f k d k 0> - y k d k l + ^f j 8 k k . (25) 
The evolution equation ( |2l] ) for T l therefore becomes 

d t r = T ik 8 j d k (3 i + p ij 8jd k k 

+0 : >d j r i ~ Vi>, f + \Pdj0i 

-2A ij 8ja + 2a(f i jk A jk + QA^d^ - ^f j 8jK). 

(26) 

In the second line we see the formula for a vector density 
of weight |, but since T l is not really a tensor density 
but is derived from Christoffel symbols we obtain extra 
terms involving second derivatives of the shift (the first 
line in the equation above). 

On the right-hand sides of the evolution equations for 
A^ and K, ([L7p and (|l8|), there occur covariant deriva- 
tives of the lapse function, and the Ricci tensor of the 
non-conformal metric. Since 

= f k l3 + 2(8 k 8jcf> + 6*184 - %7 kl 8i<t>), (27) 

where T k ij is the Christoffel symbol of the conformal met- 
ric, we have for example 

D i D i a = e-^ifldtdja - f k d k a + 2j ij 8^8 j a). (28) 

The Ricci tensor can be separated in two parts, 

R i j=R ij +Rtj, (29) 

where Rij is the Ricci tensor of the conformal metric and 
Rfj denotes additional terms depending on cf>: 

Rfj = -2A-Dj0 - 2%b k b k <t> 

+4b i <f>bj4,-4j ij b k <j>b k 4>, (30) 

with Di the covariant derivative associated with the con- 
formal metric. The conformal Ricci tensor can be written 
in terms of the conformal connection functions as 

Rij = -\*i lm 8id m % + 7 fc(i a,)f* + r k f ilj)k 

_l_~;m (^2f k l{ ifj) km + f k im f k i^j . (31) 
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A key observation here is that if one introduces the 
T l as independent variables, then the principal part of 
the right-hand side of equation ([l7]) contains the wave 
operator 7 9zd m 7ij but no other second derivatives of 
the conformal metric. This brings the evolution system 
one step closer to being hyperbolic. 

One of the reasons why we have written out the BSSN 
system in such detail is to point out a subtlety that arises 
in the actual implementation if one wants to achieve nu- 
merical stability. In the computer code we do not use the 
numerically evolved T l in all places, but follow this rule: 

• Partial derivatives 9 3 T* are computed as finite dif- 
ferences of the independent variables T l that are 
evolved using (pGj). 

• In all expressions that just require T* and not its 
derivative we substitute 7 3 ' fe r*j&(7), that is we do 
not use the independently evolved variable T l but 
recompute f 1 according to its definition (|lj) from 
the current values of 7^. 

In practice we have found that the evolutions are far less 
stable if either T l is treated as an independent variable 
everywhere, or if T* is recomputed from 7^ before each 
time step. The rule just stated helps to maintain the 
constraint T l = — 9j7 y well behaved without removing 
the advantage of reformulating the principal part of the 
Ricci tensor. 

In summary, we evolve the BSSN variables 7^, <f>, Ay, 
K, and f l according to ©, ©, @, @, and @), re- 
spectively. The Ricci tensor is separated as shown in ( p9| ) 
with each part computed according to ( |30| ) and ( |3l| ) re- 
spectively The Hamiltonian and momentum constraints 
have been used to write the equations in a particular way. 
The evolved variables T l are only used when their partial 
derivatives are needed (the one term in the conformal 
Ricci tensor ([HJ) and the advection term /3 k dkT l in the 
evolution equation for the f 1 themselves, Eq. (p6|)). 

III. THE GAUGE CONDITIONS 

We will consider families of gauge conditions that are 
not restricted to puncture data and that can be used in 
principle with any 3+1 form of the Einstein's equations 
that allows a general gauge. However, the specific fam- 
ily we test in this paper is best motivated by considering 
the BSSN system introduced above. For the present pur- 
poses, of special importance are the following two prop- 
erties of this formulation: 

• The trace of the extrinsic curvature K is treated 
as an independent variable. For a long time it has 
been known that the evolution of K is directly re- 
lated to the choice of a lapse function a. Thus, 
having K as an independent field allows one to im- 
pose slicing conditions in a much cleaner way. 



• The appearance of the "conformal connection func- 
tions" T 1 as independent quantities. As al- 
ready noted by Baumgarte and Shapiro |2l[ (see 
also the evolution equation for these quan- 

tities can be turned into an elliptic condition on 
the shift which is related to the minimal distortion 
condition. More generally, one can relate the shift 
choice to the evolution of these quantities, again al- 
lowing for a clean treatment of the shift condition. 

Our aim is to look for gauge conditions that at late 
times, once the physical system under consideration has 
settled to a final stationary state, will be able to drive 
the coordinate system to a frame where this stationar- 
ity is evident. In effect, we are looking for "symme- 
try seeking" coordinates of the type discussed by Gund- 
lach and Garfinklc and also by Brady, Creighton, and 
Thorne [^6|,0 that will be able to find the approximate 
Killing field that the system has at late times. In order 
to achieve this we believe that the natural approach is to 
relate the gauge choice to the evolution of certain com- 
binations of dynamic quantities in such a way that the 
gauge will either freeze completely the evolution of those 
quantities (typically by solving some elliptic equations), 
or will attempt to do so with some time delay (by solving 
instead parabolic or hyperbolic equations). 

We will consider the lapse and shift conditions in turn. 
Special cases of the gauge conditions that we will intro- 
duce here were recently used together with BH excision 
with remarkable results in ]l6fl , but as we will show below, 
the conditions are so powerful that in the cases tested, 
they work even without excision. 

A. Slicing conditions 

The starting point for our slicing conditions is the U K- 
freezing" condition dtK=0, which in the particular case 
when K=0 reduces to the well known "maximal slicing" 
condition. The if -freezing condition leads to the follow- 
ing elliptic equation for the lapse: 

Aa = ffdiK + n/\,,./\'< , (32) 

with A the Laplacian operator for the spatial metric 7^ . 
In the BSSN formulation, once we have solved the elliptic 
equation for the lapse, the if-freezing condition can be 
imposed at the analytic level by simply not evolving K. 

One can construct parabolic or hyperbolic slicing con- 
ditions by making either dtu or d^a proportional to dtK . 
We call such conditions "K-driver" conditions (see p8| ) . 
The hyperbolic K-driver condition has the form fi|,|18| 

d t a = ~a 2 f{a){K-K ), (33) 

where f(a) is an arbitrary positive function of a and 
Kq = K(t = 0). In our evolutions, we normally take 

/(«) = -, (34) 
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which is referred to as 1+log slicing, since empirically 
we have found that such a choice has excellent singu- 
larity avoiding properties. In Sec. fVB| we introduce a 



modification of f(a) for puncture evolutions. The hyper- 
bolic K-driver condition is in fact only a slight general- 
ization of the Bona-Masso family of slicing conditions [Q : 
d t a — ~a 2 f(a)K. 

By taking an extra time derivative of the slicing condi- 
tion above, and using the evolution equation for K, one 
can see that the lapse obeys a generalized wave equation, 



dfa = ~d t (a 2 f){K 

= a 2 /(Aa - aKijK 



K Q ) - a 2 fd t K 

lj - (3'DiK + 2af + a 2 f). (35) 



Previously we have also experimented with a somewhat 
different form of hyperbolic K-driver condition 



dfa = 



-a 2 fd t K, 



(36) 



where the right hand side vanishes in the case that K 
freezing is achieved. However, one may anticipate the 
problem that even in the case when dtK = we only 
obtain dta = const, while for (|33| ) we see that K = Kq 
implies that dta — 0. Moreover, in black hole evolutions 
where the lapse collapses to zero, condition ( |33"1) guaran- 
tees that the lapse will stop evolving, while condition d36| ) 
only implies that dta will stop evolving so the lapse can 
easily "collapse" beyond zero and become negative. For 
these reasons, in practice the condition Eq. (|33j) leads to 
more stable black hole evolutions, and this is the slicing 
condition that we consider in this paper. 

The wave speed in both cases is v a — ay/ f{a), which 
explains the need for f(a) to be positive. Notice that, 
depending on the value of f(a), this wave speed can be 
larger or smaller than the physical speed of light. This 
represents no problem, as it only indicates the speed of 
propagation of the coordinate system, i.e. it is only a 
"gauge speed". In particular, for the 1+log slicing in- 
troduced above with / = 2/ a, the gauge speed in the 
asymptotic regions (where a ~ 1) becomes v a = \/2 > 1. 
One could then argue that choosing / = 1/a should be a 
better alternative, as the asymptotic gauge speed would 
then be equal to the physical speed of light. However, 
experience has shown that such a choice is not nearly as 
robust and seems to lead easily to gauge pathologies as 
those studied in EflBOfl. 



B. Shift conditions 

In the BSSN formulation, an elliptic shift condition is 
easily obtained by imposing the "Gamma-freezing" con- 
dition d t r k =0, or using (f26|), 



2A l *d j a - 2a -f J d 3 K - Wdrf - T) k A ]h 



0. 



(37) 



Notice that, just as it happened with the if-freezing con- 
dition for the lapse, once we have solved the previous 
elliptic equations for the shift, the Gamma-freezing con- 
dition can be enforced at an analytic level by simply not 
evolving the T k . 

The Gamma-freezing condition is closely related to the 
well known minimal distortion shift condition j31J] . In 
order to see exactly how these two shift conditions are 
related, we write here the minimal distortion condition 



V,E y =0, 



(38) 



where E,,- is the so-called "distortion tensor" defined as 



(39) 



with 7y the same as before. A little algebra shows that 
the evolution equation for the conformal connection func- 



tions (|2q) can be written in terms of E^ as 



More explicitly, we have 



d t T l = 2e 



4<p 



(40) 



(41) 



We then see that the minimal distortion condition 
V J Eij = 0, and the Gamma- freezing condition dtT 1 = 
are equivalent up to terms involving first spatial deriva- 
tives of the spatial metric multiplied with the distortion 
tensor itself. In particular, all terms involving second 
derivatives of the shift are identical in both cases (but not 
so terms with first derivatives of the shift which appear 
in the distortion tensor Ey). That the difference between 
both conditions involves Christoffel symbols should not 
be surprising since the minimal distortion condition is 
covariant while the Gamma-freezing condition is not. 

Just as it is the case with the lapse, we obtain parabolic 
and hyperbolic shift prescriptions by making either dtf3 l 
or d 2 j3 l proportional to dtT 1 . We call such conditions 
"Gamma-driver" conditions. The parabolic Gamma 
driver condition has the form 



dtP 1 = F p dtT 1 



(42) 



where F p is a positive function of space and time. In 
analogy to the discussion of the hyperbolic K-driver con- 
dition there are two types of hyperbolic Gamma-driver 
conditions that we have considered, 



dfp 1 = Fd t r -ndtP 1 , 

or alternatively, 

d 2 ^ = Fd t r-( v -^)d t p\ 



(43) 



(44) 



where F and ry are positive functions of space and time. 
For the hyperbolic Gamma-driver conditions we have 
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found it crucial to add a dissipation term with coefficient 
r\ to avoid strong oscillations in the shift. Experience has 
shown that by tuning the value of this dissipation coef- 
ficient we can manage to almost freeze the evolution of 
the system at late times. 

The rational behind the two almost identical choices 
of Gamma-driver is the following. First note that if F is 
independent of time, the two choices are identical. How- 
ever, we typically choose F to be proportional to a p , 
with p some positive power (usually p = 1). Anticipating 
a collapsing lapse near the black hole this implies that 
the term FdtT 1 approaches zero and the evolution of the 
shift tends to freeze independent of the behaviour and 
numerical errors of dtT 1 . We implement the first choice 
of the Gamma-driver, ([f3|), as 

d t fi l = B\ d t B l = Fd t T l - f]B l , (45) 

and the second choice, ([5|), as 

d t f3 l = FB\ d t B l = d t f l - rjB\ (46) 

The second variant has the advantage that if F ap- 
proaches zero due to the collapse of the lapse near a 
black hole, then dtj3 l also approaches zero and the shift 
freezes. With the first variant, on the other hand, it is 
only dtB 1 = df/3 l the one that approaches zero, which 
means the shift can still evolve. Both Gamma-drivers 
can give stable black hole evolutions, although the sec- 
ond one leads to less evolution near the black holes. 

An important point that needs to be considered when 
using the hyperbolic Gamma-driver condition is that of 
the gauge speeds. Just as it happened with the lapse, 
the use of a hyperbolic equation for the shift introduces 
new "gauge speeds" associated with the propagation of 
the shift. In order to get an idea of how these gauge 
speeds behave, we will consider for a moment the shift 
condition (^3|) for small perturbations of flat space (and 
taking r)=0). From the form of dtT 1 given by equation 
( p6| ) we see that in such a limit the principal part of the 
evolution equation for the shift reduces to 

d 2 t P L = F (s jk djd k r + I S^djdk^ . (47) 

Consider now only derivatives in a given direction, say 
x. We find 

= F fo 4 + ~ S ix d x d x ^) , (48) 

which implies 

%ir = ±F%ir, (49) 

d 2 t (3i = FdlP« q^x. (50) 

We can then see that in regions where the spacetime 
is almost flat, the longitudinal part of the shift propa- 
gates with speed viong = 2yF/3 while the transverse 



part propagates with speed ^trans = VF. We therefore 
choose 

F(a) = ^a, (51) 

in order to have the longitudinal part of the shift prop- 
agate with the speed of light. The transverse part will 
propagate at a different speed, but its contribution far 
away is typically very small. 

In the next section we will turn to puncture evolutions. 
Both /(a) and F(a) will be further adjusted for the pres- 
ence of punctures. 

IV. PUNCTURES 

So far our discussion of the BSSN formulation and the 
proposed gauge conditions was quite independent of any 
particular choice of initial data, except that our gauge 
conditions are tailored for the late time stationarity of 
binary black hole mergers even though they are also ap- 
plicable in more general situations. In this Section we 
introduce puncture initial data for black holes and the 
method of puncture evolutions. 

A. Puncture initial data 

Consider the three-manifold R 3 with one or more 
points {xa, Ha, za) removed. These points we call punc- 
tures. The punctured B? arises naturally in solutions 
to the constraints in the Lichnerowicz-York conformal 
method fj2pc| ] for the construction of black hole initial 
data. In the conformal method, we introduce on the ini- 
tial hypersurface at t = the conformal variables 7y and 
Ay by 

iij = i4nu, (52) 

Aj — V'o Aiji (53) 

where ^0 is the conformal factor, and leave K untrans- 
formed. Note that Ay = ipQ 6 Aij at t = 0. 

Consider initial data with the conformally flat metric 

Hj=1p0 S i3- ( 54 ) 

Assuming that the extrinsic curvature Kij vanishes, the 
momentum constraints (m are trivially satisfied and the 
Hamiltonian constraint (Q) reduces to 

A'Vo - 0, (55) 

where A s is the flat space Laplacian. A particular solu- 
tion to this equation is 

^o = l + ^, (56) 

where r 2 = x 2 + y 2 + z 2 and we assume r ^ 0. This way 
we have obtained initial data representing a slice of a 
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Schwarzschild black hole of mass m in spatially isotropic 
coordinates on a punctured R 3 . The horizon is located at 
r = m/2. There are two asymptotically flat regions, one 
for r — ► oo and a second one at r — > 0. In fact, the metric 
is isometric under r' = m 2 /(Ar). Since ( |55| ) is linear in 
ipo one immediately obtains solutions for multiple black 
holes, for example 



BL 



1 



TOl 

2rT 



TO 2 

2»V 



(57) 



where r\ = {x — x A ) + (y - da) + (z — z A ) on an R 3 
with two punctures at (xa, Va, za), A = 1,2. These so- 
lutions were first mentioned in p3| and studied in detail 
by Brill and Lindquist in While no longer isometric, 
this initial data contains one or two black holes depend- 
ing on the separation of the punctures, but in any case 
with two separate inner asymptotic regions at the punc- 
tures. In particular, there is no physical singularity at 
the punctures, but there is a coordinate singularity at 
each puncture if one considers the unpunctured R 3 . 

Brill-Lindquist data can be generalized to longitudinal, 
non-vanishing extrinsic curvature Kij for multiple black 
holes with linear momentum and spin Q,^5|. Here one 
uses the Bowen-York extrinsic curvature, 



A*' = -^{rfPi + n>P l - (<5 y - n l n 3 ))8 k in K P l ) 
+l(n i e iU S k n l +n s e M S k ni), 



(58) 



and K = 0. The parameters P l and S l are the linear and 
angular ADM momentum, and n % — x % jr is the coordi- 
nate normal vector. The sum of two Bowen-York terms 
centered at two punctures is an explicit solution to the 
conformal momentum constraint with K = 0. 

The key observation for puncture initial data is that, 
even though there is a coordinate singularity at each 
puncture, both in the conformal factor and in the Bowen- 
York extrinsic curvature, we can rewrite the Hamiltonian 
constraint as a regular equation on R 3 without any punc- 
ture points removed. This equation possesses a unique 
solution u that is C 2 at the punctures and C°° elsewhere, 
and the original Hamiltonian constraint is solved by 



BL- 



(59) 



Working on R 3 simplifies the numerical solution of the 
constraints over methods that for example use an isomc- 
try condition at the throat of a black hole. 

Note that the puncture method for initial data can be 
applied using a conformally flat metric and the Kerr ex- 
trinsic curvature j3(|, and also to non-conformally flat 
initial data for multiple Kerr black holes |]37|-|39|]. In this 
paper we restrict ourselves to the conformally flat punc- 
ture data with Bowen-York extrinsic curvature. 



B. Puncture evolutions in the ADM system 

In this section we want to argue that one can obtain 
regular evolutions of puncture initial data without re- 
moving a region containing the puncture coordinate sin- 
gularity from the grid by, say, an isometry condition at 
the throat of the black holes as in ||, or through black 
hole excision |0,[L4j. Evolving on R 3 instead of on R 3 
with a sphere removed and an additional boundary con- 
dition imposed results in a significant technical simplifi- 
cation. 

That this is possible was first noticed experimentally 
for a single Schwarzschild puncture in [^| for ADM evolu- 
tions with singularity avoiding slicing and vanishing shift. 
By turning off the isometry condition at the throat and 
computing everywhere including next to the puncture, 
the lapse equation can still be solved for a numerically 
smooth lapse with vanishing first derivative at the punc- 
ture that collapses to zero at and around the puncture 
during the evolution. The numerical grid in these simu- 
lations is staggered around the puncture points. 

In fl, puncture evolutions are proposed as a general 
method for the evolution of the conformally flat, lon- 
gitudinal extrinsic curvature data of orb iting and spin- 
ning black holes discussed in Sec. IV A. In particular, 



an argument is given that the punctures do not evolve 
by construction. This is not a theorem about the regu- 
larity of the solutions, as is available for puncture initial 
data |34|,|35|], but it is consistent with the numerical re- 
sults. 

The basic idea is to examine the evolution equations 
and the gauge conditions at t = in the limit of small 
distance to one of the punctures. For simplicity we move 
one of the punctures onto the origin and consider the 
limit r — > 0. 

In this Section we will give a detailed version of the 
argument of 0] for the ADM equations with maximal 
slicing and vanis hing s hift, and then discuss the BSSN 
equations in Sec. IV C|. First note that for the puncture 
initial data based on (|58p and ( p9| ) we have 4>bl = 0(1 /r) 
and u — 0(1), and therefore at t = 0, 



i>o - 0(l/r), 
A 1 = 0(r), 

7y=^ = 0(l/r 4 ), 
ICy - ^ 2 A l3 = 0(l/r). 



(60) 
(61) 
(62) 
(63) 



We therefore observe that the ADM equations (0) and (^) 
for the evolution of the metric and the extrinsic curvature 
are singular at the punctures. 

The basic construction in puncture evolutions is to fac- 
tor out the time-independent conformal factor ipBL (and 
not ipo) given by the initial data, 



— tp BL Kij. 



(64) 
(65) 
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The key difference to the BSSN rescaling is that puncture 
evolutions involve a special rescaling that is independent 
of time. 

Eq. (|64| ) gives rise to a method for accurate finite dif- 
ferencing of the metric. For example, for the first partial 
derivative we have 

dklij = ^% L dkJij + %d k ip% L , (66) 

where ip% L and dki/j BL are given analytically, and 7^ and 
dkjij are assumed to remain regular during the evolution. 
By staggering the puncture between grid points one can 
therefore still obtain accurate derivatives of 7^ near the 
puncture, and this applies to all quantities derived from 
the metric and its derivatives like the Christoffcls and the 
Ricci tensor. In particular, there is no finite differencing 
across the singularity of 1/r terms. 

In general, we have tpo = u + tpBL and the analytic 
derivatives of tpo are not known, but we can still factor 
out ipBL as in ( p4| ) and ( |65| ) and obtain regular initial 
data, 

% = 1>bl<** = O(l), (67) 
K l0 =i>- B i^ 2 A l0 = 0(r 3 ). (68) 

The question is whether 7^ and K^ develop a singularity 
during the course of the evolution. 

The ADM equations for 7^ and Kij in the case of 
vanishing shift are 

8 t % = -2aK ih (69) 
d t K lj = ip^li-D^ja + aRi 3 ) 

+^ kl (k ij K kl -2K ik k jl ). (70) 

Let us examine the terms on the right hand side of the 
d t Kij equation. For 7^ = 0(1) and Kij — 0(r n ), 
the terms involving Kij are of order 0(r 2n ). Accord- 
ing to T% = f| + (r^J^-. Assuming that % 

and its derivatives are O(l), we have = O(l), but 
{T-4>bl% = 0(l/r). Hence = 0(l/r), and simi- 
larly Rij = 0(l/r 2 ). Finally, let us also assume that 
the lapse and its derivatives are of order O(l). Then 
^giDiDja = 0(r 3 ). 

With these assumptions we obtain for t = that 

d t %(0) = O(r 3 ), (71) 
d t K l3 (0) = 0(r 2 ), (72) 

where the 0(r 2 ) in the last equation is contributed by 
the term involving the Ricci tensor, the lapse terms are 
0(r 3 ), and the extrinsic curvature terms are 0(r e ). In 
order to study the time evolution, we can perform one 
finite differencing step in time from t = to t = At, for 
example, 

7« (At) - jij (0) + Atd^ (0) = 0(1), (73) 
(At) = K^ (0) + Atd t Kij (0) = 0(r 2 ). (74) 



Note that Kij has dropped from 0(r 3 ) to 0(r 2 ). How- 
ever, it is readily checked that a second finite time step 
does not further lower the order of any variable since the 
order of the right hand side in the evolution equation of 
Kij is dominated by ipg^Rij. We therefore find that 

7y(t) = 0(l), d t %{t) = 0(r 2 ), (75) 
K. l] (t) = 0{r 2 ), d t K l] {t) = 0{r 2 ). (76) 

This argument suggests that if the lapse a and its deriva- 
tives do not introduce additional singularities at the 
punctures, and if there are no singularities appearing in 
the spatial derivatives of the metric (which we have not 
completely ruled out), then the right-hand sides of the 
evolution equations for 7^ and K^ vanish at the punc- 
tures for all times. This means that 7^ and K^ should 
not evolve at all at the punctures for a regular slicing and 
vanishing shift by construction, and the same is true for 
gij and K ir 

C. Puncture evolutions in the BSSN system 

For accurate finite differencing in the BSSN system for 
puncture data we split the logarithmic conformal factor 
<f> into a singular but time-independent piece and an ad- 
ditional time-dependent contribution %, 

= X + m^BL. (77) 

It remains to be seen whether \ and the remaining 
BSSN quantities are regular throughout the evolution, 
i.e. whether the coordinate singularity can be cleanly sep- 
arated out with ipBL as in the case of ADM. To decide 
this issue we have to be specific about our gauge choice. 
In preparation for the discussion of the gauge for punc- 
ture evolutions in Sec. [v| we note some properties of the 
BSSN system near the punctures. 

Each of the BSSN variables has the following initial 
value for puncture data at t — 0, which we assume to 



evolve as indicated by the arrows: 

X = 0(1)^0(1), (78) 

K = -> 0(r 2 ), (79) 

7« = 0(1) - 0(1), (80) 

l y = 0(r 3 ) - 0(r 2 ), (81) 

r = 0^O(r). (82) 



Assume furthermore that a = O(l), and that the deriva- 
tives of the O(l) quantities are O(l). Consider now the 
following form of the evolution equations, where we have 
inserted our assumptions for a, 7y, Ay, and T l , but have 
kept the explicit dependence on (3 % , <f> and K: 

d tX = C P <P - l -aK, (83) 

d t K = CpK + 0(r 4 )(0(l) + 0{d<j>)) + \aK 2 , (84) 
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dtiij + 0(r 2 ), (85) 

dtlij = C p A l3 + KO{r 2 ) + 0{r 4 ) (0(1) 

+0(50) + 0(d 2 cp) + O(dcb) 2 ) , (86) 
d t f 1 = i),Cr' J + 0(r 2 ) + 0(r 2 )0(d(j>) 

—ardjK. (87) 

If these equations are to hold for all times, to be checked 
by time stepping as in the last section, then we require 
certain assumptions about the shift as well. In particu- 
lar, each of the terms involving Cp should be of the same 
or higher order as the other terms in the corresponding 
equation, because otherwise there could be evolution to- 
wards lower orders in r. In particular, even assuming a 
and K are regular, we have to examine the behaviour 
of Cff(j) at the puncture before we can conclude that x 
remains regular. 

Let us assume that \ an d its derivatives are regular, 
so that 

= O(lnr), d i( j> = 0(l/r), 8^ = 0(1/^). (88) 

If furthermore K — 0(r 2 ), and if the shift terms are 
of sufficiently high order, then the right hand sides of 
(p3|)-(|87|) are at least 0(r). In this case, the order of 
each equation is such that the corresponding orders in 
((78|)-(|82|) are maintained. In Sec. we show that these 
assumptions can indeed be met by a proper gauge choice, 
and hence we arrive at the statement that in this case the 
punctures do not evolve by construction. 

V. GAUGE CONDITIONS AND PUNCTURE 
EVOLUTIONS 

The main question is whether there arc lapse and shift 
conditions that behave appropriately for puncture evolu- 
tions. We will show that this is indeed the case without 
the need to introduce special boundary conditions at the 
punctures. What is required is an appropriately regular- 
ized implementation of our gauge conditions and a choice 
of initial data for lapse and shift such that the punctures 
do not evolve. 



A. Lapse for puncture evolutions 

Consider maximal slicing, which is implemented by 
choosing K — at t = and determining the lapse from 
the elliptic equation resulting from dtK = 0, which for 
vanishing shift is 

Aa = aK ij K ij . (89) 
As discussed in (7), for g. k j = ipBLliji 

Aa = ^\A^a - 8^T%d k a, (90) 



so the principal part is degenerating to zero as Q(r 4 ). 
To avoid numerical problems we therefore multiply ( |89| ) 
by V'bl' w hich normalizes the principal part but leaves 
a 0(l/r) term since = 0(l/r): 

A^a-0(l/r)d k a = 0(r 6 )a. (91) 

It turns out that standard numerical methods to solve 
this elliptic equation will find a regular solution for a for 
which dkCt vanishes sufficiently rapidly so that 0(\/r)dkU 
is zero at the puncture. Therefore, maximal slicing and 
vanishing shift lead to a sufficiently regular lapse such 
that indeed the right-hand sides of the ADM evolution 
equations for 7^ and vanish at the punctures. 

Effectively, maximal slicing implements the condition 
that the lapse has a vanishing gradient at the puncture. 
Notice that this condition is very different to an isometry- 
type condition, where the lapse would be forced to be — 1 
at the puncture. Fig. [l] shows a schematic representation 
in the Kruskal diagram of the type of slices one obtains 
in the case of a single Schwarzschild BH when using three 
different boundary conditions for the lapse (while keeping 
the same interior slicing condition): odd at the throat, 
even at the throat and zero gradient at the puncture. 
When looking at these plots it is important to keep in 
mind that the puncture corresponds to a compactifica- 
tion of the second asymptotically flat region, and is at 
an infinite distance to the left of the plots. Notice that, 
in all three cases, far away on the right hand side of the 
plot the slices approach the Schwarzschild slices (in fact, 
if we use maximal slicing and ask for the lapse to be odd 
we recover the Schwarzschild slices everywhere). Also, in 
the case with an odd lapse the slices do not penetrate the 
horizon, but in the other two cases they do. 

Since maximal slicing is computationally expensive, we 
often use 1+log slicing that mimics the behaviour of the 
maximal lapse in that it also is singularity avoiding and 
the lapse drops to zero when the physical singularity ap- 
proaches. Analytically, however, the 1+log lapse does 
not necessarily drop to zero at the puncture. Starting 
with a = 1 and K = everywhere at t = 0, we see from 
the evolution equations for the lapse and for K , 

d t a^~a 2 f(a)(K-K ) 1 (92) 

d t K = f3 i d i K +\aK 2 + 0(r 3 ), (93) 

that neither a nor K evolve at the puncture if the shift 
or the derivative of K vanishes at the puncture. That 
means that the lapse will remain 1 at the puncture, and 
the inner asymptotically flat region will evolve. Numer- 
ically, one may expect this to be problematic if there is 
not sufficient resolution, as will be normally the case. The 
code can become unstable near the puncture, and even 
if it remains stable the event horizon may not prevent 
numerical inaccuracies to evolve into the outer regions. 
While this may happen, in practice the 1+log lapse does 
collapse, apparently precisely because of a lack of resolu- 
tion, and the code remains stable. Even in this case we 
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Odd lapse Even lapse 




evolution in the still poorly resolved region between the 
puncture and the horizon, which is why we do not use 
this option routinely. Instead of guessing an initial lapse 
that minimizes the amount of initial evolution one should 
use the lapse (and shift) derived from a quasi-equilibrium 
thin-sandwich puncture initial data set, which however is 
currently not available. 



Zero gradient at puncture 




FIG. 1. Schematic representation on the Kruskal diagram 
of the effect of the different boundary conditions on the slices 
obtained. The first panel shows the case of an odd lapse at 
the throat, the second panel the case of an even lapse at the 
throat, and the last panel the case of a lapse with zero gradient 
at the puncture. The dashed lines show the singularities and 
the dotted lines the event horizon. 



typically obtain convergence in the outer regions where 
the lapse did not collapse. 

In this paper we experiment with 1+log slicing with 
f(a) = 2/a replaced by f(a,ipBL) = ^^bl/ -^ so that 



d t a 



-2a^ L (K-K ), 



(94) 



i.e. we have introduced a factor that for m > can drive 
the lapse to zero at the puncture. For m = we ob- 
tain standard 1+log slicing. A natural choice is m = 4 
since then the singularity in i/j^l exactly matches the de- 
generacy of the principal part of the second order wave 
equation associated with the lapse evolution, see ( |35| ) and 
(|90|). In particular, for m — 4 the wave speed is regular 
at the puncture. 

Both choices of 1+log slicing with initial lapse equal to 
unity have been found to lead to stable evolutions of black 
holes with a lapse that satisfies the regularity conditions 
assumed in the previous sections. Another approach to 
obtain a vanishing lapse at the puncture is to start with 
a different initial lapse, for example 



a(t = 0) = i' B 2 L = 0(r 2 ), 



(95) 



so that the lapse is zero at the puncture initially and 
there is no evolution due to a non- vanishing lapse at the 
puncture. The power —2 is chosen so that the initial lapse 
has the same limit for r — > oo as the lapse anisotropic = 
(1 — §?)/(! + 5r) °f the static Schwarzschild metric in 
isotropic coordinates. In practice, we have found that 
with such an initial lapse there sometimes is too much 



B. Shift for puncture evolutions 

For long term stable evolutions, we want to construct 
a shift condition that counters slice-stretching. However, 
for arbitrary non- vanishing shift, equations (83-8?]) show 
that the punctures will evolve. It is possible to have 
the punctures move across the grid because of a non- 
vanishing shift. One problem would be the numerical 
treatment of the coordinate singularity at the punctures, 
which so far was based on analytic derivatives of the time- 
independent conformal factor ipBL- While a solution of 
this problem may be possible, we focus here on finding a 
shift condition that counters slice-stretching while simul- 
taneously satisfying a fall-off condition for the shift such 
that the punctures do not evolve at all when using the 
BSSN equations. 

As a first step it is instructive to consider (3 % — rri 1 — 
x l = 0(r) (with n % a radial unit vector) near the punc- 
ture. In this case several terms in the Lie derivatives 
cancel exactly and we have 



C rn K = x k d k K. 



(96) 
(97) 
(98) 



However, 



C rn <j> = x k d kX + x k 8 k ln^si + ~ = 0(1), (99) 

so x wn l evolve without a special combination of x, K , 
and a. Furthermore, the Lie derivative will not be as 
simple if the shift is not exactly spherically symmetric. 
We therefore turn to 



ft = 0(r 3 ). 



(100) 



This happens to be the necessary condition (assuming 
integer powers of r) for the norm of the shift in the non- 
conformal metric to be zero at the puncture 



7^/3^' = 0(l/r 4 )«JyW 



With ft 



r 2 x l we now have 



£r»n% = 0{r 2 ){x k d k % + %) = 0(r 2 ), 
C r s n (j) = 0(r 2 ), 



(101) 

(102) 
(103) 



All other Lie derivative terms also turn out to be of or- 
der 0{r 2 ). Finally, the shift derivatives in the evolution 
equation for T l are 
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BjCtf* = d,0(r 2 ) = 0(r). (104) 

In this sense the evolution of T l poses the strictest con- 
dition on the fall-off of the shift. 

The question remains how we guarantee the 0(r 3 ) fall- 
off in the actual shift condition. We can enforce such a 
fall-off by choosing 

p i (t = 0) = 0, (105) 

and by changing the coefficient F in the hyperbolic 
Gamma-driver to 

%^L) = |«feL = 0(r"), (106) 

where we typically choose n = 2 or 4. Note that the 
two versions of the Gamma-driver differ by the term 
d t F/F = d t a/a, which in the case of 1+log slicing equals 
—2ip' l g L (K — K a ). Let us ignore the diffusion term. If the 
shift has evolved into (3 l = 0(r 3 ), then d t f l = 0(r). 
With n > 2 we have 

d 2 t (3 l = 0(r 3 ). (107) 

In fact, we have found that in the case of evolutions of just 
one black hole (distorted or not), changing F is not really 
necessary since there is enough symmetry in the problem 
to guarantee that the shift will have the correct fall-off at 
the puncture even without introducing the extra factor of 
iPbl into F. When dealing with two black holes, however, 
this is no longer the case and the factor ips^L is required 
to stop the shift from evolving at the punctures. 

Anticipating the numerical results, let us point out 
that due to a lack of numerical resolution the shift of- 
ten looks like 0(r) even though we have chosen n = 2 or 
4. Somewhat surprisingly, even in these cases the gauge 
conditions are able to approximately freeze the evolution, 
for which we do not yet have a good analytic explanation. 

C. Vanishing of the shift at the punctures 

Combining our choice of puncture initial data with the 
lapse and shift conditions above, we have the expectation 
that the BSSN variables will not evolve at the punctures. 
This is a natural situation considering that the punctures 
represent an asymptotically flat infinity, and that there 
is no linear momentum at the inner infinities. Perform- 
ing the transformation r — > 1/r for puncture data with 
the Bowen-York extrinsic curvature defined in d5S| ) shows 
that the 1/r 3 spin terms are mapped to 1/r 3 terms, but 
the 1/r 2 linear momentum terms are mapped to 1/r 4 
terms and there are no 1/r 2 terms, and therefore there 
is no linear momentum at the inner infinity. In other 
words, viewed from the other asymptotic end, the black 
hole does not move in the data we use. 

One can add a 1/r 4 term to to make the holes move 
in the inner ends, but then the puncture initial data con- 
struction and the puncture evolutions on R 3 fail for a 



lack of regularity at the punctures. In general, with a 
different choice of extrinsic curvature that does not sat- 
isfy the fall-off conditions of the Bowen-York data J5^), 
there can be non-trivial or even singular evolution at the 
punctures in both the ADM and the BSSN systems. 

In summary, our puncture initial data corresponds to 
two black holes which are momentarily at rest at their in- 
ner asymptotic ends. For a given coordinate system the 
black holes could start moving if there is a non- vanishing 
shift at the punctures, but we explicitly construct a van- 
ishing shift at the punctures. The main consequence for 
puncture data of orbiting and colliding black holes is that 
by construction the inner asymptotic ends of the black 
hole will not move in our coordinate system, i.e. the punc- 
tures remain glued to the grid. That still allows for gen- 
eral dynamics around the punctures, which shows in the 
evolution of the metric and extrinsic curvature. For ex- 
ample, the apparent horizon can easily grow, drift and 
change shape, but it can not cross over the punctures 
for geometrical reasons, since the apparent horizon area 
would become infinite before they could do that. For or- 
biting black holes, since the punctures do not move by 
construction, it seems natural to combine the method of 
puncture evolutions with a corotating coordinate system 
to minimize the evolution of the metric data. We leave 
this option for future work. 

VI. NUMERICS 

The numerical time integration in our code uses an 
iterative Crank-Nicholson scheme with 3 iterations, see 
e.g. [Q. Derivatives are represented by second order 
finite differences on a Cartesian grid. We use standard 
centered difference stencils for all terms, except in the ad- 
vection terms involving the shift vector (terms involving 
^di). For these terms we use second order upwind along 
the shift direction. We have found the use of an upwind 
scheme in such advection-type terms crucial for the sta- 
bility of our code. Notice that this is the only place in our 
implementation where any information about causality is 
used (i.e. the direction of the tilt in the light cones). 

A. Outer boundary condition 

At the outer boundary we use a radiation (Sommer- 
feld) boundary condition. We start from the assump- 
tion that near the boundary all fields behave as spherical 
waves, namely we impose the condition 

f = fo + u(r-vt)/r. (108) 

Where /o is the asymptotic value of a given dynamical 
variable (typically 1 for the lapse and diagonal metric 
components, and zero for everything else), and v is some 
wave speed. If our boundary is sufficiently far away one 
can safely assume that the speed of light is 1, so v = 1 
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for most fields. However, the gauge variables can eas- 
ily propagate with a different speed implying a different 
value of v. 



In practice, we do not use the boundary condition (10S 
as it stands, but rather we use it in differential form: 



d t f + vd r f-v(f-f o )/r = 0. 



(109) 



Since our code is written in Cartesian coordinates, we 
transform the last condition to 



■d t f + vdif+^(f-f o ) = 0. 



(110) 



We finite difference this condition consistently to second 
order in both space and time and apply it to all dynamic 
variables (with possible different values of /o and v) at 
all boundaries. 

There is a final subtlety in our boundary treatment. 
Wave propagation is not the only reason why fields evolve 
near a boundary. Simple infall of the coordinate ob- 
servers will cause some small evolution as well, and such 
evolution is poorly modeled by a propagating wave. This 
is particularly important at early times, when the radia- 
tive boundary condition introduces a bad transient ef- 
fect. In order to minimize the error at our boundaries 
introduced by such non-wavelike evolution, we allow for 
boundary behavior of the form: 



f = fo + u(r-vt)/r + h(t)/r n , 



(111) 



with h a function of t alone and n some unknown power. 
This leads to the differential equation 

dtf + vd r f - - (/ - fo) = (1 - nv) + 



h'(t) 



for large r, (112) 



or in Cartesian coordinates 



^d t f + vd if+ V -^(f-f )^^. (113) 

This expression still contains the unknown function 
h'it). Having chosen a value of n, one can evaluate the 
above expression one point away from the boundary to 
solve for h'(t), and then use this value at the boundary 
itself. Empirically, we have found that taking n = 3 al- 
most completely eliminates the bad transient caused by 
the radiative boundary condition on its own. 



stretch limited resources as far as possible, is to intro- 
duce a radial coordinate transformation that decreases 
the resolution with distance. Such coordinate transfor- 
mations can also be applied to 3D Cartesian grids, see 
the "fish-eye transformation" in P,pT| . 

In order to make the outer boundary conditions as sim- 
ple as possible, we would like for the resolution to be 
constant at the location of the outer boundaries. That 
is, we want constant high resolution in the region con- 
taining the black holes, then we want a region where the 
resolution decreases with distance and finally we want 
a region (containing the outer boundaries) with constant 
low resolution. Denoting the physical radius by r and the 
coordinate radius by r c , the previous requirements can be 
met with the following radial coordinate transformation 



r = ar c + (1 — a) 
— In I cosh 



2tanh(^) 
r c ~ r 



In cosh ■ 



ro 



(114) 



where a is a parameter specifying the scale factor in grid 
spacing, ro is the radius of the transition region and s is 
the width of the transition region. 

By differentiating r in equation (114) with respect to 
r c we find that dr/dr c = 1 for r c = and dr/dr c = a 
for r c — > oo as required. Note that the radial r coordi- 
nate is mapped to ar c plus a non- vanishing constant, and 
therefore the Jacobian of this transformation does not 
correspond to just a simple rescaling of radial resolution. 
The transformation (114) we refer to as the "transition 
fish-eye" . 

An important point to keep in mind when using a fish- 
eye transformation is the fact that both the asymptotic 
values of metric components and the physical speed of 
light (and gauge speeds) will be affected by the transfor- 
mation. This means that special care should be taken 
when applying boundary conditions. 



VII. APPLICATIONS 

In the numerical application of our method we focus 
on establishing the basic validity of the puncture evolu- 
tions with the hyperbolic shift. We consider evolutions 
of the spherically symmetric Schwarzschild spacetime, of 
a single distorted black hole, and of the head-on colli- 
sion of two Brill-Lindquist punctures. We will report on 
orbiting binary systems elsewhere. 



B. Fish-eye transformation 

Setting up a reasonable numerical simulation, there is 
always the conflicting interest of having the boundary as 
far out as possible and having as good resolution as pos- 
sible. With limited numerical resources it is almost never 
possible to obtain both at the same time. One way to 



A. Evolution of a single Schwarzschild puncture 

For the evolution of Schwarzschild we use the Cartoon 
method of ]4l| ] for implementing axisymmetric systems 
with 3D Cartesian finite differencing stencils. Choosing 
the z-axis as the axis of symmetry, we evolve a 3D Carte- 
sian slab with just 5 points in the y-direction. On the 
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FIG. 2. Schwarzschild black hole evolved for t = 1000M. 
Shown are lapse a and shift component j3 x along the x-axis, 
which are (anti-)symmetric about x — 0. By that time lapse 
and shift are approximately static. The lapse has collapsed to 
zero at the puncture and approaches one in the outer region. 
The shift crosses zero at the puncture, pointing away from 
the puncture and thereby halting the infall of points towards 
the puncture. 



y = plane standard 3D stencils are computed, and the 
data at the points with y ^ are obtained by interpola- 
tion in the x-direction in the y — plane and by tensor 
rotation about the z-axis. For Schwarzschild wc also use 
the reflection symmetry in the z = plane. 

We choose the Schwarzschild puncture data of 



Sec. [V A with m = 1.0M and the apparent horizon at 
r = 0.5M. As we have discussed, there are several choices 
for the gauge conditions. For the Schwarzschild puncture, 
we initialize lapse and shift to a — 1 and (3 l — 0. We con- 
sider 1+log slicing, (33), and the hyperbolic shift, (44), 
with the specific choice of 



2a- 1 1>% L , F=-a^ 2 L , V = 2.0/M. (115) 



In Fig. U we show lapse and shift for an evolution with 
201 points in x- and z-direction, starting at the staggered 
point at the origin and extending to about 20M with a 
grid spacing of 0.1M. We plot the data after an evolution 
of t = 1000A/, which corresponds to 40000 time steps 
with a Courant factor of 0.25. 

Lapse and shift show the characteristic feature of punc- 
ture evolutions. Both lapse and shift are zero at the 
puncture, indicating that there is no evolution at the in- 
ner asymptotically flat end of the black hole. The lapse 
approaches one in the outer region, while the shift points 
outward from the puncture and approaches zero in the 
outer region. The shift counters the infall of points to- 
wards the puncture, thereby stopping the slice stretching. 

Fig. U shows various other quantities of the same 




5 10 

X(M) 

FIG. 3. Schwarzschild black hole evolved for t = 1000M. 
Shown are the BSSN variables (f>, K, y xx , A xx , and T x along 
the x-axis, and also the Hamiltonian constraint H . 



Schwarzschild run at t = 1000M. Note the behaviour 
near the puncture, which at this resolution appears to be 
regular but is not sufficiently resolved. 

In Fig. |we compare data from this run with a run for 
identical parameters except that instead of Eq. (Q) we 
use Eq. @ with j) = 2.8/ M for the shift. The differences 
are quite small in the case of this Schwarzschild run. 

In Fig. H we show the maximum of the shift and the 
root-mean-square value of the Hamiltonian constraint as 
a function of time. After a short time interval of less 
than 100 M (recall that previous runs with vanishing 
shift lasted only to about 30-40M!), evolution is approx- 
imately frozen for more than 3000M. The observed drift 
in various quantities is crucially affected by the value of rj 
that determines the diffusion in the hyperbolic Gamma- 
driver. In Fig. |B| we compare again the two versions of 
the Gamma-driver, and note that two different values of 
77, 2.0/M and 2.8/M, are used to obtain long term stabil- 
ity. It is a matter of experimentation to find a suitable 
value of r\ in dependence on the various parameters in 
the run. Runs may crash before 100M for a bad choice 
of 77. On the other hand, once determined for a particu- 
lar initial data set and set of grid parameters, we found 
that the runs were rather robust under small variations. 
It would be useful to have a dynamic determination and 
adaptation of ?y, but this is currently not available. 

Having established the basic features and the stabil- 
ity of our gauge choice, we want to study convergence 
for Schwarzschild. A crucial issue is whether we ob- 
tain convergence near the puncture. We choose three 
grid sizes and resolutions: 201 points in both the x- and 
z-directions for a resolution of 0.025M, 401 points at 
0.0125M, and 801 points at 0.00625M. With a Courant 
factor of 0.25 in the BSSN evolution scheme it takes 160, 
320, and 640 iterations, respectively, for an evolution 
time of 1M. The outer boundary is at about 5M. We 
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FIG. 4. Schwarzschild black hole evolved for t = 1000M. 
Shown is a comparison along the x-axis between two versions 
of the hyperbolic Gamma-driver for the shift, Eq. (^3j) (dashed 
line) and Eq. ( p4| ) (solid line). 
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FIG. 5. Schwarzschild black 

hole evolved for t = 3000 M. Shown are the maximum of 
the shift and the root-mean-square value of the Hamiltonian 
constraint as a function of time, again for two versions of the 
hyperbolic Gamma-driver for the shift, Eq. (^) (dashed line) 
and Eq. ( p4| ) (solid line), with diffusion parameter r\ — 2.0/ 'M 
and 7] = 2.8/ M, respectively. After a short time interval dur- 
ing which lapse and shift adjust themselves dynamically, the 
evolution slows down significantly. 
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FIG. 6. Schwarzschild black hole evolved to t — 5M at 
three high resolutions, demonstrating second order conver- 
gence at the puncture. 



choose the same gauge as in (115), except that in F we 
use ipgL instead of for a broader profile of the shift 
near the puncture. 

Fig. ^| shows the Hamiltonian constraint along the x- 
axis near the single Schwarzschild puncture at the three 
resolutions, rescaled by the corresponding factors ex- 
pected for second order convergence. The coincidence of 
the three lines indicates clean second order convergence. 
Therefore, for such high resolutions the BSSN system 
exhibits the expected regular and convergent behaviour 
near the punctures. It is remarkable that even at a four 
times coarser resolution of 0.1M the evolutions remain 
stable. 

Note in particular that the shift in Fig. || seems to 
be linear at the puncture, in contrast with the expected 
0(r 3 ) behaviour. Fig. |?] shows the effect of different pow- 
ers of ipBL hr the shift equation for the grid parameters of 
the medium resolution run of the convergence test. We 
use the shift equation (0), and 



F=-a^l, V = 2.0/M, 



(116) 



with different values for n. Fig. shows the shift for 
t = 1M. Note the resolution that is required to make 
the 0(r 3 ) behaviour visible for n > 2. By t = 10M, 
the shift for n — 2 is no longer completely resolved at 
the puncture with a grid spacing of 0.0125M, but as we 
have seen, even at coarser resolutions the approximate 
0(r) behaviour of the shift at the puncture allows stable 
evolutions. 



B. Evolution of a single, distorted black hole 

The second application we present is that of a distorted 
BH. Referring to f42]] , we choose a distortion parameter 
Qo = 0.5, position 770 = 0, and width a = 1. The ADM 
mass of this system is M = 1.83. Such data has been 
previously evolved in 2D and in 3D using excision. Here 
we discuss a 3D puncture evolution with octant symme- 
try, 129 3 points and a grid spacing of 0.1M — 0.183. The 
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FIG. 7. Schwarzschild black hole evolved for t = 1M. 
Shown is the effect of varying the power n in Vb" m ^ ne 
shift equation for /3 X along the x-axis. 

outer boundary is at about 12. 8M. For the gauge we use 
1+log slicing with the initial lapse not unity but given 
by (^5|) and the hyperbolic shift condition (^4|) with 

/ = 2a~ Vk. F = \a^ B 2 Ll r, = 1.25/Af » 0.68. (117) 

In Fig. ||, we show the evolution of the lapse and the 
shift component (3 X along the x-axis. Note that the shift, 
although vanishing initially, develops the needed profile 
simply through its evolution equation, without any spe- 
cial initial condition. After a short while, the evolution 
effectively freezes, allowing the waves to propagate on an 
effectively fixed BH background, just as one would like. 

For comparison, we show in Fig. ^ the evolution of 
the radial component of the metric, Jrr/ipBL' f° r the 
new gauge condition (lower panel) and for a singularity 
avoiding slicing run with 1+log slicing and vanishing shift 
(upper panel). For 1+log slicing and vanishing shift we 
see the well-known slice stretching effect. With the new 
gauge evolution is slowed significantly at late times. The 
peak of the metric near x = 0.5M grows to about 12 by 
time t = 20M and does not grow significantly after that 
until t — 400M (lower panel), while for vanishing shift 
already at t = 30M the peak in the metric has reached 
40 without any sign of slowing growth (upper panel). 

For the new gauge we expect that we can reliably ex- 
tract the waveform for the ring-down, and this is indeed 
the case as shown in Fig. [l0| 

C. Head-on collision of two Brill-Lindquist punctures 

The third application we present is that of a head 
on collision of two Brill-Lindquist BH's. The param- 
eters for these simulations are mi = 777,2 = 0.5Af, 
Ci = {0,l.r515M, 0}, C 2 = {0,-1.1515M,0}, where 
777i and 77i2 are the masses of the BH's and Ci and C 2 
are the locations of the two punctures. These parame- 
ters correspond to an initial separation of the BH's equal 
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FIG. 8. Lapse and shift for the evolution of a single dis- 
torted BH. After around 20M, the evolution of lapse and shift 
slows down significantly (note the time labels). The approach 
to the final profile in lapse and shift is not monotonic. 

to that of an approximate ISCO configuration as deter- 
mined in p3| . Such data has been previously evolved 
without shift with the Lazarus technique that combines 
short, fully numerical evolutions with a close limit ap- 
proximation for the wave extraction Q (see |l2] ] for runs 
starting at larger separation). 

We present two types of runs for the head-on collision 
starting at the approximate ISCO separation. In the first 
type we use 1+log slicing and the hyperbolic Gamma- 
driver ( ^3|) with 

/-2a- 1 , F=^st »7 = 2.8/M. (118) 

with an initial lapse equal to one and an initial shift equal 
to zero. We also use the transition fish-eye with param- 
eters a = 3, s = 1.2M and tq = 5.5M. This places the 
outer boundary at a distance of 25. 8M with central res- 
olutions 0.128M, 0.064M and 0.032M and gridsizes 96 3 , 
192 3 and 384 3 respectively in octant mode. 

In Fig. [ll] we show the extracted 1 = 2 and m = 
waveforms until t = 80M for all three resolutions. The 
code actually continued beyond t = 140M at the highest 
resolution (more than t = 200M at the lower resolutions) 
before we stopped it, due to the fact that it was computa- 
tionally fairly expensive. Initially there seem to be some 
small amplitude oscillations superposed on the larger os- 
cillations. These seem to be related to an initial wave 
pulse in the lapse moving outwards as the lapse collapses 
from its initial value, which is not quite handled by the 
wave extraction procedure. However these oscillations 
decrease with increasing resolution. With / = 2a~ 1 tp% L 
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FIG. 9. The evolution of the radial metric function 
r )rr/ip%L f° r a distorted BH along the a;— axis. The upper 
panel shows the slice stretching in the metric for singularity 
avoiding slicing with vanishing shift, while the lower panel 
shows the metric for the new gauge conditions. Without shift 
the metric grows out of control after t — 40M, while with 
the new shift condition a peak begins to form initially but 
later almost freezes as lapse and shift drive the BH into an 
essentially static configuration (note the time labels). 



as we used in the previous examples instead of ( |118[ ) , the 
oscillations are larger, probably because the lapse is more 
dynamic in the initial p hase of the evolution. But as al- 
ready mentioned in Sec. V A, even with / = 2a~ 1 ^% L the 
lapse collapses at the punctures. After about t = 80M 
we see some non-quasi-normal features in the waveform, 
that are most probably due to contaminations from the 
outer boundary. 

For a wave signal A extracted at three resolutions, A, 
2 A and 4A the order of convergence, a, can be estimated 
as 



log 2 



A(4A) - A(2A) 



A(2A) - A(A) 



(119) 



In Fig. [12|we show this estimate of the convergence fac- 
tor for the three waveforms from Fig. [0]. Several features 
in this figure deserve comment. First of all, for the first 
15M the signal is very small, so the estimate of the con- 
vergence order is not very accurate. Secondly, the phase 
evolution of the waveforms is somewhat resolution depen- 
dent. This means that the curves cross over each other at 
different times, leading to the spikes clearly visible in the 
plot. The differences in phase evolution seem to decrease 
with increasing resolution, although only at somewhere 
between first and second order. However excluding the 
initial part and the spikes, we see a reasonable second 
order convergence in the waveforms up to t = 80M. 
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The solid line shows the 
extracted at a radius of 5.45M for the even-parity distorted 
BH described in the text, while the two dashed lines show 
the result of the same simulation carried out in the 2D and 
3D code with vanishing shift. The 2D code crashes at around 
t = 100M and the 3D code crashes around t = 40M. The 
lower panel shows a fit for the time interval from t = 9A/ 
to t = 80M to the two lowest quasi normal modes of the 
BH for the new gauge conditions in 3D, confirming that the 
ring-down of the distorted BH is simulated accurately. 

In Table |, we try to circumvent the problem of the 
differently evolving phase by locating the extrema of the 
waveforms and estimating the convergence order using 
these extremal values even if they do not occur at the 
same time. As can be seen, except for the first maximum, 
there is generally nice second order convergence in the 
amplitude. In the case of the first maximum, it can be 
seen from Fig. |ll| that the difference between the three 
resolutions is very small and that especially the lowest 
resolution is influenced by the pulse in the lapse moving 
out. 

As a second type of gauge choice we use maximal slic- 
ing and the hyperbolic gamma driver condition with the 



Extremum 
1 
2 
3 
4 
5 
6 



log 2 |(A(4A) - A(2A))/(A(2A) - A(A))\ 
1.17 
2.11 
2.00 
1.95 
1.96 
2.24 



TABLE I. The convergence of the amplitude for the first 
six local extrema of the extracted £ — 2, m — waveforms for 
the head on collision of two Brill-Lindquist BH's extracted at 
radius 18. 5M. 
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FIG. 11. The extracted £ = 2, m = waveforms for the 
head on collision of two Brill-Lindquist BH's at three different 
resolutions extracted at 18.5M. The resolutions for the solid, 
dashed and dash-dotted line are 0.032M, 0.064M and 0.128M 
respectively. 
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FIG. 12. The convergence order in the extracted 1 = 2, 
m = waveforms for the head on collision of two 
Brill-Lindquist BH's extracted at 18.5M, based on the same 
three resolutions (0.032M, 0.064M, 0.128M) as in Fig. [n[ 

same shift parameters as in the 1+log case, except for 
the fact that rj = 2.0/M here. In this case the resolution 
is 0.128M and the grid size is 80 3 in octant mode. The 
fish-eye parameters are a = 4, s = 1.2M and ro = 5.0M 
placing the outer boundary at a distance of 26M. This 
run ran for about a month on two processors on a dual 
1.7 GHz Xeon workstation reaching t = 5224M, until the 
machine went down due to an unrelated problem. By 
that time, the evolution was almost completely frozen as 
can be seen from Fig. [l3| showing the common apparent 
horizon mass as function of time. Most of the evolution 
occurs before t = 200M and after that there is just a slow 
drift of the apparent horizon mass giving about 10% error 
at t = 5000A/. 

In Fig. [l4| we plot the extracted waveform with a loga- 
rithmic time scale (actually ln(f + 1)) in order to be able 
to see the features in the beginning of the waveform, while 
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FIG. 13. The apparent horizon mass for the head on colli- 
sion of two Brill-Lindquist BH's with maximal slicing. 
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FIG. 14. The extracted t = 2, m = waveforms for the 
head on collision of two Brill-Lindquist BH's with maximal 
slicing. Note that we plot t + 1 in order to be able to use a 
logarithmic scale on the time axis. 

still showing that it is constant and very close to zero at 
t = 5000A/ . The features in the initial part of this wave- 
form are very similar to the features in the 1+log run 
of the same resolution. However it is completely smooth 
in the initial phase, where the 1+log waveform has the 
small amplitude oscillations, since with an elliptic lapse 
condition, there is no wave pulse in the lapse moving 
outwards. 

As mentioned before, these evolutions were done in oc- 
tant symmetry. We repeated the maximal slicing evolu- 
tion in bitant symmetry (reflection about one coordinate 
plane) , with exactly the same physical and gauge param- 
eters. However, this evolution died at about t — 280Af, 
showing some clearly asymmetric features in the lapse 
and metric components in the directions where the sym- 
metry is not imposed. We first encountered such a depen- 
dence of the stability of the BSSN system on the octant 
symmetry in excision runs of a single black hole pq ] . The 
current result supports the conclusion that the stability 
problem is not directly linked to the excision technique or 
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the gauge conditions, but is probably intrinsic to BSSN. 
We are currently investigating the cause of this problem. 

In conclusion, with the new gauge conditions we can 
evolve not only single black hole systems but also the 
head on collision of two black holes with dynamically 
adjusting lapse and shift and reach an almost static so- 
lution for the final black hole. While we have argued in 
detail why the punctures should not evolve, and while it 
is plausible that there is sufficient freedom in the gauge 
to almost freeze the evolution of a single, spherical black 
hole, it is remarkable that the method is successful even 
in the region close to and between two black holes. 

VIII. CONCLUSION 

We have discussed a new family of coordinate con- 
ditions for 3D numerical relativity that are powerful, 
efficient, easy to implement, and respond naturally to 
the spacetime dynamics. An application of these condi- 
tions to previously difficult BH spacetimes shows their 
strength: even without excision, they allow distorted 
and colliding BH spacetimes to be evolved for more than 
two orders of magnitude longer than possible previously, 
for thousands of M rather than tens of M, while keep- 
ing errors down to a few percent and allowing accurate 
waveform extraction. The evolution methods and gauge 
choices discussed here have already passed preliminary 
tests for orbiting punctures. Work is in progress to mod- 
ify the shift condition for corotating coordinates. 
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